## This script will re-run the models used to impute baseline temperature and PM data

library(HEIfunctions)
library(spBayes)
library(splines)
library(xtable)
library(pscl)
library(DMwR)
set.seed(51)

dir = ''

## Load data with missing temperature (if you want to impute it here)
load('datforimpute.RData')
dat = datforimpute

##############################################################################
########### IMPUTE MISSING TEMPERATURE DATA (If not already done) #################
##############################################################################
indices <- which(is.na(dat$temperature_annual_max_mean) == T)
imputetemp <- knnImputation(cbind(dat$temperature_annual_max_mean, dat$Longitude, dat$Latitude), 3)
dat$temperature <- dat$temperature_annual_max_mean
dat$temperature[indices] <- imputetemp[indices,1]
dat$temporig = dat$temperature_annual_max_mean  ## this line added on 9/4/2015 when recreating to prevent a later (insignificant) error
save(dat, file = 'datforimpute_withtemp.Rda')
set.seed(51)
#############################################################
############## BASELINE SPATIAL PREDICTIONS #################
#############################################################

#### CHOOSE WHICH BASELINE VARIABLE TO IMPUTE ######
dat$pm10base = dat$pm10base1990

obsdat 	<- dat[which(is.na(dat$pm10base)==F),]
misdat 	<- dat[which(is.na(dat$pm10base)==T),]

varprior 	<- c(.15,1)
samp <- sqrt(rigamma(10000,varprior[1], varprior[2]))

obscoords 	<- cbind(obsdat$Longitude, obsdat$Latitude)
miscoords 	<- cbind(misdat$Longitude, misdat$Latitude)
phiest 	<- makephi(obscoords, scale=10) #1.12531
phiprior 	<- c(makephi(obscoords, scale=3), makephi(obscoords, scale=50)) #(.3376, 5.6265)


starting 	<- list("phi"=phiest,"sigma.sq"=10, "tau.sq"=10)
tuning 	<- list("phi"=0.01, "sigma.sq"=0.1, "tau.sq"=0.1)
priors 	<- list("phi.Unif"=phiprior, "sigma.sq.IG"=varprior,"tau.sq.IG"=varprior)

nsamp		<- 50000
nburn		<- 10000

formula = pm10base~ a  + Median_income + HS_rate + Urban_rate + White_cens + log(POPULATION) + temperature

summary(lm(formula, data=obsdat))

outmodel <- spatpred(formula, data=obsdat, data.p = misdat,  
                     coords = obscoords, coords.p = miscoords, nburn=nburn,
                     starting=starting, tuning=tuning, priors = priors,
                     cov.model="exponential", n.samples = nsamp,
                     verbose = T)

save(outmodel, file = "20150427pm1990preds.Rda")

